Primero cargamos los paquetes que vamos a utilizar
library(tidyverse)
library(knitr)
library(here)
También los paquetes específicos para EDA, pero esta vez no lo mostramos en el documento final porque en el chunk ponemos echo = FALSE:
Y alguno más para modelos:
library(broom)
library(ggfortify)
library(GGally)
library(modelsummary) #- remotes::install_github('vincentarelbundock/modelsummary')
library(equatiomatic) #- remotes::install_github("datalorax/equatiomatic")
library(cluster)
library(factoextra)
library(FactoMineR)
library(rpart)
library(rattle)
Hemos visto en tutoriales anteriores como cargar, arreglar y manejar datos con R à la tidyverse. En este tutorial vamos a ver algunos paquetes y funciones que nos pueden ayudar a acelerar el análisis exploratorio inicial de nuestros datos.
R es un lenguaje/entorno que fue creado para hacer estadística y casi no hemos visto nada de modelos estadísticos. Al final de este tutorial veremos como estimar algunos modelos con R. Tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo (solo algo) de Machine Learning.
Una de las máximas de un data scientist es, como señalan aquí, “Know Your Data”. Además nos explican que:
Data exploration is the art of looking at your data, rapidly generating hypotheses, quickly testing them, then repeating again and again and again. The goal of data exploration is to generate many promising leads that you can later explore in more depth.
With exploratory data analysis, you’ll combine visualisation and transformation with your curiosity and scepticism to ask and answer interesting questions about data
Simplemente dos recursos:
En este repo de Github, Mateusz Staniak mantiene un listado de paquetes y recursos relacionados con la EDA.
Por último, un libro, ay de 2016, pero de Roger Peng: Exploratory Data Analysis with R, covers visualization, graphics devices, dimension reduction, and data wrangling.
Estaría bien que vieseis alguno de lo screencast de David Robinson.
Son datos de bebes valencianos para los años 2010 y 2017. Los datos provienen del INE.Concretamente están aquí.
Cargamos los datos y el diccionario en memoria de R (Global Environment)
df_orig <- rio::import(here::here("datos", "nacidos", "df_partos_23v_CV_10_17.rds")) #- 90.236 x 23
diccionario <- rio::import(here::here("datos", "nacidos", "dicc_partos_23v_CV_10_17.xlsx"))
Los datos originales tienen 23 variables de 90236 partos.
Vamos a restringir el análisis a partos de un sólo bebe
df <- df_orig %>% filter(MULTSEX %in% c(1,2)) #- solo partos con 1 bebe (88.176)
Con lo que ahora tenemos 88176 registros
Además voy a quitar algunas variables:
df <- df %>% select(MUNPAR.n, SEMANAS, CESAREA.f, PESON1, SEXO1.f, NUMHV, EDADM.2, EDADP.2, ESTUDIOM.f, ESTUDIOP.f, PAISNACM.n, PAISNACP.n, fecha_parto, fecha_parto_anterior) #- me quedo con 14 variables
Finalmente, nuestro df tiene 14 variables de 88176 partos.
En este apartado haremos repaso de algunas funciones y paquetes que nos pueden ser de utilidad en el análisis exploratorio inicial.
Conviene saber acceder a los nombres de las variables. Por ejemplo con base::names()
names(df)
#> [1] "MUNPAR.n" "SEMANAS" "CESAREA.f"
#> [4] "PESON1" "SEXO1.f" "NUMHV"
#> [7] "EDADM.2" "EDADP.2" "ESTUDIOM.f"
#> [10] "ESTUDIOP.f" "PAISNACM.n" "PAISNACP.n"
#> [13] "fecha_parto" "fecha_parto_anterior"
Y saber cambiarlo si lo necesitamos
names(df)[2] <- "nuevo_nombre_variable_2"
names(df)
#> [1] "MUNPAR.n" "nuevo_nombre_variable_2"
#> [3] "CESAREA.f" "PESON1"
#> [5] "SEXO1.f" "NUMHV"
#> [7] "EDADM.2" "EDADP.2"
#> [9] "ESTUDIOM.f" "ESTUDIOP.f"
#> [11] "PAISNACM.n" "PAISNACP.n"
#> [13] "fecha_parto" "fecha_parto_anterior"
Ejercicio: por favor, volved a dejar el nombre de la segunda variable como estaba: SEMANAS
names(df)[2] <- "SEMANAS"
names(df)
#> [1] "MUNPAR.n" "SEMANAS" "CESAREA.f"
#> [4] "PESON1" "SEXO1.f" "NUMHV"
#> [7] "EDADM.2" "EDADP.2" "ESTUDIOM.f"
#> [10] "ESTUDIOP.f" "PAISNACM.n" "PAISNACP.n"
#> [13] "fecha_parto" "fecha_parto_anterior"
Si necesitamos arreglar los nombres, la función janitor::clean_names() es fantástica para esto. Lo hace automáticamente, además tiene opciones para hacerlo como más te guste.
df_jugar <- df
janitor::clean_names(df_jugar) %>% names()
#> [1] "munpar_n" "semanas" "cesarea_f"
#> [4] "peson1" "sexo1_f" "numhv"
#> [7] "edadm_2" "edadp_2" "estudiom_f"
#> [10] "estudiop_f" "paisnacm_n" "paisnacp_n"
#> [13] "fecha_parto" "fecha_parto_anterior"
Siempre-siempre hay que saber de que tipo son nuestras variables. Esto se puede hacer de muchas formas. Os muestro algunas:
str(df) #- muestra la estructura interna de un objeto R
#> Classes 'tbl_df', 'tbl' and 'data.frame': 88176 obs. of 14 variables:
#> $ MUNPAR.n : chr "Villajoyosa/Vila Joiosa, la" "Alzira" "Gandia" "Valencia" ...
#> $ SEMANAS : num 39 40 41 34 40 38 40 40 38 37 ...
#> $ CESAREA.f : Factor w/ 2 levels "Sin cesárea",..: 1 2 1 2 1 1 1 1 1 1 ...
#> $ PESON1 : num 3000 3260 3650 1720 2800 ...
#> $ SEXO1.f : Factor w/ 2 levels "Varón","Mujer": 1 1 1 2 1 1 1 2 2 1 ...
#> $ NUMHV : num 0 0 0 2 2 0 0 0 0 0 ...
#> $ EDADM.2 : num 34.4 21.3 24.6 35.2 40.9 ...
#> $ EDADP.2 : num NA NA NA NA NA NA NA NA NA NA ...
#> $ ESTUDIOM.f : Factor w/ 3 levels "Primarios","Medios",..: NA NA NA NA NA NA NA NA NA NA ...
#> $ ESTUDIOP.f : Factor w/ 3 levels "Primarios","Medios",..: NA NA NA NA NA NA NA NA NA NA ...
#> $ PAISNACM.n : chr "Bulgaria" "Bulgaria" "Bulgaria" "Bulgaria" ...
#> $ PAISNACP.n : chr NA NA NA NA ...
#> $ fecha_parto : POSIXct, format: "2017-05-01" "2017-07-01" ...
#> $ fecha_parto_anterior: POSIXct, format: NA NA ...
inspectdf::inspect_types(df) #- muestra de que tipo son las variables
#> # A tibble: 4 x 4
#> type cnt pcnt col_name
#> <chr> <int> <dbl> <list>
#> 1 numeric 5 35.7 <chr [5]>
#> 2 factor 4 28.6 <chr [4]>
#> 3 character 3 21.4 <chr [3]>
#> 4 POSIXct POSIXt 2 14.3 <chr [2]>
Yo suelo utilizar estas dos funciones
df_aa <- pjp.funs::my_df_estadisticos_basicos(df) #- estadisticos basicos del df
df_bb <- pjp.funs::my_df_valores_unicos(df) #- valores únicos de cada variable de df
La función summarytools::dfSummary() da un resumen muy útil
#library(summarytools)
view(dfSummary(df)) #- genera un fichero con un resumen muy útil y agradable de ver para cada variable
Otra alternativa es usar skimr::skim(df)
skimr::skim(df)
Siempre hay que ver la existencia de NA’s en nuestro df. Tres paquetes que nos ayudan en esta tarea son DataExplorer, visdat y naniar.
DataExplorer::plot_missing(df)
naniar::gg_miss_var(df, show_pct = TRUE) #- nº de % de NA's en cada columna
El paquete visdat es útil para ver los NA’s pero el siguiente chunk no nos va a funcionar. A ver si conseguís que funcione:
visdat::vis_dat(df)
Como resulta que visdat se queja de que hay muchas variables, entonces, vamos a visualizar sólo las variables que sí tienen NA’s:
Ejercicio: seleccionar las variables que tienen algún NA.
zz_con_NAs <- df %>% select_if(anyNA)
zz_factores <- df %>% select_if(is.factor)
zz_x <- df %>% select(PESON1, EDADM.2, EDADP.2, SEXO1.f)
Visualicemos las variables con NA's
visdat::vis_dat(zz_con_NAs)
naniar::gg_miss_upset(df)
naniar: Concretamente veremos la ocurrencia de NA’s para los grupos/categorías de la variable ESTUDIOM.f1naniar::gg_miss_var(df, facet = ESTUDIOM.f, show_pct = TRUE) #- faceted por la variable ESTUDIOM.f
rm(list = ls(pattern = "^zz")) #- borra los objetos cuyo nombre empiece por zz
# rm(zz_con_NAs, zz_factores, zz_x)
Pues no está claro, depende … pero vamos a dejar un df SIN NAs; y vamos a aprovechar para practicar algo de manejo de datos con R:
Ejercicio: quita los NA’s del df
zz_1 <- df[complete.cases(df), ] #- con base-R
zz_2 <- df %>% filter(complete.cases(.)) #- con dplyr y base-R
zz_3 <- df %>% tidyr::drop_na() #- con tidyr
Ok, si quitamos todos los registros/partos que tienen algún NA en alguna de las 14 variables nos quedamos con 26492 registros.
PERO, la última variable, fecha_parto_anterior tiene muchos NA’s porque para muchas mujeres ese parto/bebe es el primero, así que vamos a refinar la eliminación de NA’s
zz_4 <- df %>% tidyr::drop_na(-fecha_parto_anterior) #- quito filas con NA sin tener en cuenta fecha_parto_anterior
zz_5 <- df %>% tidyr::drop_na(1:13 ) #- quito filas con NA en las 13 primeras variables
Sí, vamos a hacerlo:
df <- df %>% tidyr::drop_na(-fecha_parto_anterior)
rm(list = ls(pattern = "^zz")) #- borra los objetos cuyo nombre empiece por zz
Ahora nuestro df tiene 57051 partos, registros o filas.
A veces es útil esta función
visdat::vis_expect(df, ~.x == 1) #- ver la prevalencia de un determinado valor (puede ser de utilidad, mira la ayuda!!)
valores_a_buscar <- c("Con cesárea", "Universidad")
visdat::vis_expect(df, ~.x %in% valores_a_buscar)
Ejercicio: seleccionar los partos en los que ambos padres tienen estudios universitarios y el parto se realizó con cesárea.
zz <- df %>% filter(ESTUDIOM.f == "Universidad" &
ESTUDIOP.f == "Universidad" &
CESAREA.f == "Con cesárea" )
Vamos a analizar un poco las variables numéricas.
Ejercicio: crea un df solo con las variables numéricas
df_numeric <- df %>% select_if(is.numeric)
Por ejemplo, summarytools::descr() nos ayuda a calcular rápidamente estadísticos descriptivos de las variables numéricas.
summarytools::descr(df, style = "rmarkdown")
N: 57051
| EDADM.2 | EDADP.2 | NUMHV | PESON1 | SEMANAS | |
|---|---|---|---|---|---|
| Mean | 32.30 | 34.98 | 0.62 | 3237.90 | 38.96 |
| Std.Dev | 5.27 | 5.80 | 0.80 | 508.86 | 1.89 |
| Min | 13.01 | 14.59 | 0.00 | 450.00 | 20.00 |
| Q1 | 29.19 | 31.44 | 0.00 | 2960.00 | 38.00 |
| Median | 32.69 | 34.86 | 0.00 | 3250.00 | 39.00 |
| Q3 | 35.94 | 38.36 | 1.00 | 3560.00 | 40.00 |
| Max | 56.04 | 71.46 | 11.00 | 6350.00 | 45.00 |
| MAD | 4.95 | 5.19 | 0.00 | 444.78 | 1.48 |
| IQR | 6.75 | 6.93 | 1.00 | 600.00 | 2.00 |
| CV | 0.16 | 0.17 | 1.29 | 0.16 | 0.05 |
| Skewness | -0.38 | 0.25 | 2.06 | -0.61 | -2.31 |
| SE.Skewness | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 |
| Kurtosis | 0.15 | 1.11 | 9.60 | 2.35 | 10.78 |
| N.Valid | 57051.00 | 57051.00 | 57051.00 | 57051.00 | 57051.00 |
| Pct.Valid | 100.00 | 100.00 | 100.00 | 100.00 | 100.00 |
zz <- summarytools::descr(df)
Y DataExplorer lo vamos a utilizar para obtener histogramas o gráficos de densidad para las variables numéricas
DataExplorer::plot_histogram(df)
DataExplorer::plot_density(df)
Veamos la matriz de correlaciones entre las variables numéricas:
df %>% select_if(is.numeric) %>% corrr::correlate() %>% gt::gt()
| SEMANAS | PESON1 | NUMHV | EDADM.2 | EDADP.2 | |
|---|---|---|---|---|---|
| SEMANAS | NA | 0.52092275 | -0.05170781 | -0.04590809 | -0.02587726 |
| PESON1 | 0.52092275 | NA | 0.06702148 | 0.00018377 | 0.02165019 |
| NUMHV | -0.05170781 | 0.06702148 | NA | 0.23533870 | 0.23102964 |
| EDADM.2 | -0.04590809 | 0.00018377 | 0.23533870 | NA | 0.64879388 |
| EDADP.2 | -0.02587726 | 0.02165019 | 0.23102964 | 0.64879388 | NA |
Otra forma de calcular y visualizar las correlaciones entre las variables de un df
DataExplorer::plot_correlation(df)
DataExplorer::plot_correlation(df, type = 'c') #- solo las continuas
Otra forma de visualizar las correlaciones:
zz <- df %>% inspectdf::inspect_cor()
df %>% inspectdf::inspect_cor() %>% inspectdf::show_plot()
Para ver diferencias en la distribución de las variables numéricas frente a una v. categórica es muy útil DataExplorer::plot_boxplot(). Por ejemplo frente a la variable ESTUDIOM.f:
DataExplorer::plot_boxplot(df, by = "ESTUDIOM.f")
O, por ejemplo frente a la variable CESAREA.f
my_vv <- names(df)[3]
my_vv <- "CESAREA.f"
DataExplorer::plot_boxplot(df, by = my_vv)
df %>% explore::explore(EDADM.2, ESTUDIOM.f, target = CESAREA.f)
df %>% select(EDADM.2, PESON1, ESTUDIOM.f, CESAREA.f) %>% explore::explore_all(target = CESAREA.f)
DataExplorer::plot_scatterplot() nos permite hacer rápidamente un scatterplot de todas las variables del df frente a una variable, por ejemplo el peso del bebe: PESON1
my_vv <- "PESON1"
DataExplorer::plot_scatterplot(df, by = my_vv, sampled_rows = 1000L)
Si necesitásemos un df con todas las variables categóricas, ¿cómo lo obtendrías?
df_categoricas <- df %>% select_if(negate(is.numeric)) #- !!!!
Para visualizar rápidamente los valores de las variables categóricas, tenemos inspectdf::inspect_cat(), que nos devuelve un gráfico con la distribución de todas las variables categóricas.
inspectdf::inspect_cat(df) %>% inspectdf::show_plot(high_cardinality = 1)
DataExplorer::plot_bar()DataExplorer::plot_bar(df)
GGally: una extensión de ggplot2 que sirve para muchas cosas, por ejemplo:df_jugar <- df %>% select(EDADM.2, SEXO1.f, SEMANAS, PESON1)
GGally::ggpairs(df_jugar, progress = FALSE)
GGally::ggscatmat(df_jugar, columns = c("EDADM.2", "SEMANAS", "PESON1"))
burro:burro attempts to make EDA accessible to a larger audience by exposing datasets as a simple Shiny App
#- devtools::install_github("laderast/burro")
burro::explore_data(df, outcome_var = colnames(df))
explore.Instead of learning a lot of R syntax before you can explore data, the explore package enables you to have instant success. You can start with just one function - explore() - and learn other R syntax later step by step
explore(df)
df %>% explore(CESAREA.f)
df %>% explore(EDADM.2, ESTUDIOM.f, target = CESAREA.f)
zz <- df %>% janitor::get_dupes(PESON1) #- janitor::get_dupes() calcula los valores de PESON1 repetidos
zz <- df %>% janitor::remove_empty(c("rows", "cols")) #- quita variables y filas vacias
dplyr::distinct() es muy útil para ver las filas con valores únicos. Por ejemplo:zz <- df %>% distinct(ESTUDIOM.f, CESAREA.f)
zz %>% kable()
| ESTUDIOM.f | CESAREA.f |
|---|---|
| Primarios | Con cesárea |
| Universidad | Sin cesárea |
| Primarios | Sin cesárea |
| Medios | Sin cesárea |
| Universidad | Con cesárea |
| Medios | Con cesárea |
zz %>% kable(format = "html") %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
| ESTUDIOM.f | CESAREA.f |
|---|---|
| Primarios | Con cesárea |
| Universidad | Sin cesárea |
| Primarios | Sin cesárea |
| Medios | Sin cesárea |
| Universidad | Con cesárea |
| Medios | Con cesárea |
df %>% count(PAISNACM.n) %>%
mutate(PAISNACM.n = fct_reorder(PAISNACM.n, n)) %>%
arrange(desc(n)) %>% slice(1:10) %>%
ggplot(aes(PAISNACM.n, n)) + geom_col() + coord_flip() +
labs(title = "Nº de bebes por nacionalidad de la madre")
Muchas veces hay que presentar resultados básicos como una tabla de casos o de frecuencias. para ello tenemos muchas posibilidades:
summarytoolsEl paquete summarytools.
summarytoolsis an R package providing tools to neatly and quickly summarize data. It can also make R a little easier to learn and to use, especially for data cleaning and preliminary analysis.
freq() provee tablas con conteos y frecuenciassummarytools::freq(df$SEXO1.f, style = "rmarkdown")
Type: Factor
| Freq | % Valid | % Valid Cum. | % Total | % Total Cum. | |
|---|---|---|---|---|---|
| Varón | 29534 | 51.77 | 51.77 | 51.77 | 51.77 |
| Mujer | 27517 | 48.23 | 100.00 | 48.23 | 100.00 |
| <NA> | 0 | 0.00 | 100.00 | ||
| Total | 57051 | 100.00 | 100.00 | 100.00 | 100.00 |
summarytools::ctable(df$SEXO1.f, df$CESAREA.f)
Cross-Tabulation, Row Proportions
SEXO1.f * CESAREA.f
Data Frame: df
| CESAREA.f | Sin cesárea | Con cesárea | Total | |
| SEXO1.f | ||||
| Varón | 20472 (69.3%) | 9062 (30.7%) | 29534 (100.0%) | |
| Mujer | 19486 (70.8%) | 8031 (29.2%) | 27517 (100.0%) | |
| Total | 39958 (70.0%) | 17093 (30.0%) | 57051 (100.0%) |
Incluso se puede hacer un test chi-cuadrado
summarytools::ctable(df$ESTUDIOM.f, df$CESAREA.f, chisq = TRUE)
Cross-Tabulation, Row Proportions
ESTUDIOM.f * CESAREA.f
Data Frame: df
| CESAREA.f | Sin cesárea | Con cesárea | Total | |
| ESTUDIOM.f | ||||
| Primarios | 15734 (73.2%) | 5769 (26.8%) | 21503 (100.0%) | |
| Medios | 11764 (68.8%) | 5335 (31.2%) | 17099 (100.0%) | |
| Universidad | 12460 (67.5%) | 5989 (32.5%) | 18449 (100.0%) | |
| Total | 39958 (70.0%) | 17093 (30.0%) | 57051 (100.0%) |
| Chi.squared | df | p.value |
|---|---|---|
| 168.1 | 2 | 0 |
names(df)
[1] “MUNPAR.n” “SEMANAS” “CESAREA.f”
[4] “PESON1” “SEXO1.f” “NUMHV”
[7] “EDADM.2” “EDADP.2” “ESTUDIOM.f”
[10] “ESTUDIOP.f” “PAISNACM.n” “PAISNACP.n”
[13] “fecha_parto” “fecha_parto_anterior”
janitorLa verdad es que janitor es un paquete fantástico!!!
janitor has simple functions for examining and cleaning dirty data. It was built with beginning and intermediate R users in mind and is optimized for user-friendliness. Advanced R users can already do everything covered here, but with janitor they can do it faster and save their thinking for the fun stuff.
df %>% janitor::tabyl(CESAREA.f) %>% gt::gt()
| CESAREA.f | n | percent |
|---|---|---|
| Sin cesárea | 39958 | 0.7003909 |
| Con cesárea | 17093 | 0.2996091 |
df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f) %>% gt::gt()
| CESAREA.f | Primarios | Medios | Universidad |
|---|---|---|---|
| Sin cesárea | 15734 | 11764 | 12460 |
| Con cesárea | 5769 | 5335 | 5989 |
df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f, SEXO1.f) #%>% gt::gt()
#> $Varón
#> CESAREA.f Primarios Medios Universidad
#> Sin cesárea 7969 6004 6499
#> Con cesárea 3102 2752 3208
#>
#> $Mujer
#> CESAREA.f Primarios Medios Universidad
#> Sin cesárea 7765 5760 5961
#> Con cesárea 2667 2583 2781
Las tablas con R-base tampoco están mal. El problema que tienen es que los resultados no se almacenan en data.frames:
table(df$CESAREA.f, df$ESTUDIOM.f, useNA = "always")
#>
#> Primarios Medios Universidad <NA>
#> Sin cesárea 15734 11764 12460 0
#> Con cesárea 5769 5335 5989 0
#> <NA> 0 0 0 0
my_tabla <- table(df$CESAREA.f, df$ESTUDIOM.f)
my_tabla %>% as.data.frame() %>% gt::gt()
| Var1 | Var2 | Freq |
|---|---|---|
| Sin cesárea | Primarios | 15734 |
| Con cesárea | Primarios | 5769 |
| Sin cesárea | Medios | 11764 |
| Con cesárea | Medios | 5335 |
| Sin cesárea | Universidad | 12460 |
| Con cesárea | Universidad | 5989 |
zz <- prop.table(my_tabla,1)
mosaicplot(t(zz), color = TRUE, main = "% de Cesáreas para niveles educativos de la madre")
Evidentemente R permite implementar múltiples técnicas y modelos estadísticos, así como hacer múltiples y variados contrastes de hipótesis. En esta pagina web tenéis un buena introducción a estos tema. Veamos algún ejemplo: por ejemplo:
t.test(df$PESON1)
t.test(df$PESON1 ~ df$CESAREA.f)
t.test(df$PESON1 ~ df$SEXO1.f) %>% psycho::analyze()
t.test(df$PESON1 ~ df$SEXO1.f, var.equal = TRUE, conf.level = .95) %>% psycho::analyze() %>% summary()
cor_results <- cor.test(df$PESON1, df$EDADM.2) # Compute a correlation and store its result
psycho::analyze(cor_results) # Run the analyze function on the correlation
#> The Pearson's product-moment correlation between df$PESON1 and df$EDADM.2 is not significant, very small, and positive (r(57049.00) = 0.00, 95% CI [-0.01, 0.01], p > .1).
library(purrr)
library(broom)
df %>% group_by(ESTUDIOM.f) %>%
summarise(t_test = list(t.test(PESON1))) %>%
mutate(tidied = map(t_test, tidy)) %>%
tidyr::unnest(tidied) %>%
ggplot(aes(estimate, ESTUDIOM.f)) + geom_point() +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
labs(y = "")
chi-squared test: https://data-flair.training/blogs/chi-square-test-in-r/chisq.test(df$CESAREA.f, df$SEXO1.f)
#>
#> Pearson's Chi-squared test with Yates' continuity correction
#>
#> data: df$CESAREA.f and df$SEXO1.f
#> X-squared = 15.155, df = 1, p-value = 0.00009901
chisq.test(df$CESAREA.f, df$ESTUDIOM.f)
#>
#> Pearson's Chi-squared test
#>
#> data: df$CESAREA.f and df$ESTUDIOM.f
#> X-squared = 168.07, df = 2, p-value < 0.00000000000000022
En este post, Jonas Kristoffer Lindeløv, nos presenta un cuadro con los contrastes de hipótesis más frecuentes y cómo efectuar esos contrastes en el contexto de los modelos de regresión con R.
En este apartado vamos a ver como estimar modelos lineales y no lineales con R. Para una introducción a la estimación de modelos en R puedes ir aquí.
Tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo de Machine Learning.
Antes restrinjamos un poco el df:
df_m <- df %>% select(PESON1, SEMANAS, SEXO1.f, CESAREA.f, EDADM.2, EDADP.2, ESTUDIOM.f, ESTUDIOP.f)
La función para estimar modelos lineales es lm(), así que vamos a utilizarla para estimar nuestro primer modelo con R. Estimaremos un modelo lineal con variable a explicar el peso del bebe (PESON1) en función de todas las demás variables en df_m.
mod_1 <- lm(PESON1 ~ . , data = df_m)
summary(mod_1)
Generalmente las variables categóricas se introducen en los modelos mediante variables dummies. Crear dummies en R es sencillo, solo tienes que tener los datos como factor o como texto y R creará las dummies por ti cuando introduzcas la variable en lm(). Eso sí, es más facil elegir la categoría de referencia si la variable es un factor.
Para entender como fija los regresores para las variables categóricas, mira esto:
levels(df$SEXO1.f)
#> [1] "Varón" "Mujer"
levels(df$ESTUDIOM.f)
#> [1] "Primarios" "Medios" "Universidad"
Si necesitas cambiar la categoría de referencia siempre puedes usar forcast::fct_relevel()
zz <- forcats::fct_relevel(df$SEXO1.f, "Mujer")
levels(zz)
#> [1] "Mujer" "Varón"
Veamos que hay en el objeto mod_1
str(mod_1)
#listviewer::jsonedit(mod_2, mode = "view") ## Interactive option
Lógicamente, a veces querremos seleccionar las variables explicativas:
mod_1 <- lm(PESON1 ~ SEMANAS, data = df)
mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1.f + CESAREA.f + EDADM.2 , data = df_m)
mod_1 <- lm(log(PESON1) ~ log(SEMANAS), data = df_m)
summary(mod_1)
#>
#> Call:
#> lm(formula = log(PESON1) ~ log(SEMANAS), data = df_m)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.80542 -0.08174 0.00610 0.08778 1.39651
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.63860 0.04359 14.65 <0.0000000000000002 ***
#> log(SEMANAS) 2.02921 0.01190 170.47 <0.0000000000000002 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.1476 on 57049 degrees of freedom
#> Multiple R-squared: 0.3375, Adjusted R-squared: 0.3375
#> F-statistic: 2.906e+04 on 1 and 57049 DF, p-value: < 0.00000000000000022
Si queremos introducir interacciones entre los regresores, podemos usar el operador :, aunque casi mejor hacerlo directamente con I():
mod_1 <- lm(PESON1 ~ SEMANAS + SEMANAS:EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*EDADM.2), data = df_m)
summary(mod_1)
#>
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * EDADM.2), data = df_m)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2633.6 -279.6 -14.5 262.9 3709.3
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -2252.786335 37.656076 -59.83 < 0.0000000000000002 ***
#> SEMANAS 138.961398 0.994026 139.80 < 0.0000000000000002 ***
#> I(SEMANAS * EDADM.2) 0.060934 0.008869 6.87 0.00000000000647 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 434.2 on 57048 degrees of freedom
#> Multiple R-squared: 0.272, Adjusted R-squared: 0.2719
#> F-statistic: 1.066e+04 on 2 and 57048 DF, p-value: < 0.00000000000000022
Si queremos introducir las variables originales y también las interacciones entre ellas, podemos hacerlo directamente o o utilizar el operador *:
mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2 + SEMANAS:EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*ESTUDIOM.f, data = df_m)
summary(mod_1)
#>
#> Call:
#> lm(formula = PESON1 ~ SEMANAS * ESTUDIOM.f, data = df_m)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2630.8 -280.8 -15.4 263.9 3727.1
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -2287.0303 57.9094 -39.493 <0.0000000000000002
#> SEMANAS 141.5391 1.4842 95.366 <0.0000000000000002
#> ESTUDIOM.fMedios -29.0793 90.4969 -0.321 0.7480
#> ESTUDIOM.fUniversidad 183.7079 91.1839 2.015 0.0439
#> SEMANAS:ESTUDIOM.fMedios 0.9445 2.3193 0.407 0.6838
#> SEMANAS:ESTUDIOM.fUniversidad -4.0701 2.3383 -1.741 0.0818
#>
#> (Intercept) ***
#> SEMANAS ***
#> ESTUDIOM.fMedios
#> ESTUDIOM.fUniversidad *
#> SEMANAS:ESTUDIOM.fMedios
#> SEMANAS:ESTUDIOM.fUniversidad .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 434.2 on 57045 degrees of freedom
#> Multiple R-squared: 0.2719, Adjusted R-squared: 0.2718
#> F-statistic: 4260 on 5 and 57045 DF, p-value: < 0.00000000000000022
Recuerda que si queremos introducir algún regresor que sea la multiplicación de dos variables, tendremos que hacerlo con I()
mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*SEMANAS), data = df_m)
summary(mod_1)
#>
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * SEMANAS), data = df_m)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2635.6 -282.7 -15.6 261.5 3803.8
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -3116.4624 255.2162 -12.211 < 0.0000000000000002 ***
#> SEMANAS 188.5366 13.8452 13.617 < 0.0000000000000002 ***
#> I(SEMANAS * SEMANAS) -0.6514 0.1878 -3.469 0.000522 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 434.3 on 57048 degrees of freedom
#> Multiple R-squared: 0.2715, Adjusted R-squared: 0.2715
#> F-statistic: 1.063e+04 on 2 and 57048 DF, p-value: < 0.00000000000000022
Cambien puede sernos de utilidad la función poly()
lm(PESON1 ~ poly(SEMANAS, degree=3), data = df_m)
#>
#> Call:
#> lm(formula = PESON1 ~ poly(SEMANAS, degree = 3), data = df_m)
#>
#> Coefficients:
#> (Intercept) poly(SEMANAS, degree = 3)1
#> 3238 63314
#> poly(SEMANAS, degree = 3)2 poly(SEMANAS, degree = 3)3
#> -1507 -7131
Una vez que sabemos la sintaxis de stats::lm(), vamos a ver como podemos acceder a la información que devuelve lm(). ya hemos visto que los resultados se almacenan en una lista, así que podemos utilizar los operadores habituales de las listas para acceder a la información. Por ejemplo:
zz_betas <- mod_1[[1]]
zz_residuals <- mod_1[[2]]
zz_residuals <- mod_1[["residuals"]]
zz_residuals <- mod_1$residuals
zz_betas_2 <- mod_1[1] #- veis la diferencia
zz_betas_3 <- mod_1[[1]] #- veis la diferencia
zz_betas_2a <- mod_1["coefficients"]
zz_betas_3a <- mod_1$"coefficients"
Afortunadamente ya hay construidas algunas funciones para manipular/ver los resultados de estimación. Por ejemplo:
library(ggfortify)
summary(mod_1) #- tabla resumen
summary(mod_1)$coefficients #- tabla resumen con los coeficientes
coefficients(mod_1) #- coeficientes estimados
confint(mod_1, level = 0.95) #- Intervalos de confianza para los coeficientes
#fitted(mod_1) #- predicciones (y-sombrero, y-hat)
#residuals(mod_1) #- vector de residuos
#model.matrix(mod_1) #- extract the model matrix
anova(mod_1) #- ANOVA
vcov(mod_1) #- matriz de var-cov para los coeficientes
#influence(mod_1) #- regression diagnostics
# diagnostic plots
layout(matrix(c(1,2,3,4),2,2)) #- optional 4 graphs/page
plot(mod_1) #-
library(ggfortify)
autoplot(mod_1, which = 1:6, ncol = 2, colour = "steelblue")
Para hacer predicciones con observaciones fuera de la muestra has de usar la función predict, proporcionándole las observaciones a predecir como un df.
nuevas_observaciones <- df_m %>% slice(c(3,44,444))
predict(mod_1, newdata = nuevas_observaciones) #- predicciones puntuales
predict(mod_1, newdata = nuevas_observaciones, type = 'response', se.fit = TRUE) #- tb errores estándar predictions
predict(mod_1, newdata = nuevas_observaciones, interval = "confidence") #- intervalo (para el valor esperado)
predict(mod_1, newdata = nuevas_observaciones, interval = "prediction") #- intervalo (para valores individuales)
Hay muchas más funciones para valorar la idoneidad de un modelo. Por ejemplo aquí.
El paquete GGally permite hacer:
GGally::ggcoef(mod_1)
df_jugar <- df_m %>% select(EDADM.2, SEXO1.f, SEMANAS)
GGally::ggpairs(df_jugar)
Lo habitual es que al estimar un modelo, tengamos situaciones de heterocedasticidad, clustering etc … Podemos ajustar los errores estándar a estas situaciones. El paquete de referencia solía ser sandwich, aunque quizás ahora el paquete de referencia sea estimatr. Por ejemplo se pueden obtener errores estándar robustos con estimatr::lm_robust()2.
#- install.packages("emmeans")
mod_1 <- lm(PESON1 ~ SEMANAS, data = df_m)
mod_1_ee <- estimatr::lm_robust(PESON1 ~ SEMANAS, data = df_m)
broom::tidy(mod_1_ee, mod_1, conf.int = T)
#> term estimate std.error statistic
#> 1 (Intercept) -2240.6797 60.139009 -37.25834
#> 2 SEMANAS 140.6182 1.533625 91.69008
#> p.value
#> 1 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003194318
#> 2 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
#> conf.low conf.high df outcome
#> 1 -2358.5525 -2122.8069 57049 PESON1
#> 2 137.6123 143.6242 57049 PESON1
broomUn opción interesante es utilizar el paquete broom. Este paquete tiene 3 funciones útiles:
tidy()
augment()
glance()
zz <- broom::tidy(mod_1, conf.int = TRUE)
zz %>% gt::gt()
| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | -2240.6797 | 37.6300489 | -59.54496 | 0 | -2314.4348 | -2166.9246 |
| SEMANAS | 140.6182 | 0.9647186 | 145.76088 | 0 | 138.7274 | 142.5091 |
mod_1 %>% broom::glance() %>% select(adj.r.squared, p.value)
#> # A tibble: 1 x 2
#> adj.r.squared p.value
#> <dbl> <dbl>
#> 1 0.271 0
broom::glance(mod_1)
#> # A tibble: 1 x 11
#> r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
#> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
#> 1 0.271 0.271 434. 21246. 0 2 -4.27e5 8.55e5 8.55e5
#> # ... with 2 more variables: deviance <dbl>, df.residual <int>
zz <- broom::augment(mod_1)
broommod_1 %>% tidy() %>% filter(p.value < 0.05)
#> # A tibble: 2 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) -2241. 37.6 -59.5 0
#> 2 SEMANAS 141. 0.965 146. 0
broom, pero antes vamos a recordar alguna cosa de ggplot2:ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1, color = ESTUDIOM.f)) +
geom_point(alpha = 0.1) + geom_smooth(method = "lm")
ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1, color = SEXO1.f)) +
geom_point(alpha = 0.1) + geom_smooth(method = "lm")
ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1, color = SEXO1.f)) +
geom_point(alpha = 0.1) + geom_smooth()
Ahora sí viene el ejemplo en que se usa el paquete broom
mod_1 <- lm(PESON1 ~ EDADM.2 + SEXO1.f , data = df_m)
summary(mod_1)
#>
#> Call:
#> lm(formula = PESON1 ~ EDADM.2 + SEXO1.f, data = df_m)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2841.04 -280.59 19.01 318.31 3058.89
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 3288.96818 13.29270 247.427 <0.0000000000000002 ***
#> EDADM.2 0.07408 0.40164 0.184 0.854
#> SEXO1.fMujer -110.83978 4.23833 -26.152 <0.0000000000000002 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 505.8 on 57048 degrees of freedom
#> Multiple R-squared: 0.01185, Adjusted R-squared: 0.01181
#> F-statistic: 342 on 2 and 57048 DF, p-value: < 0.00000000000000022
td_mod_1 <- mod_1 %>% broom::augment(data = df_m)
td_mod_1 %>% ggplot(mapping = aes(x = EDADM.2, y = PESON1, color = SEXO1.f)) +
geom_point(alpha = 0.1) +
geom_line(aes(y = .fitted, group = SEXO1.f))
td_mod_1 %>% ggplot(mapping = aes(x = EDADM.2, y = .fitted, color = SEXO1.f)) +
geom_line(aes(group = SEXO1.f))
(!!!!) Por ejemplo también permite fácilmente estimar modelos por grupos:
df_m %>% group_by(SEXO1.f) %>% do(tidy(lm(PESON1 ~ SEMANAS, .)))
#> # A tibble: 4 x 6
#> # Groups: SEXO1.f [2]
#> SEXO1.f term estimate std.error statistic p.value
#> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 Varón (Intercept) -2363. 52.1 -45.4 0
#> 2 Varón SEMANAS 145. 1.34 109. 0
#> 3 Mujer (Intercept) -2158. 53.2 -40.5 0
#> 4 Mujer SEMANAS 137. 1.36 100. 0
O hacer gráficos de los intervalos de los coeficientes:
td_mod_1 <- tidy(mod_1, conf.int = TRUE)
ggplot(td_mod_1, aes(estimate, term, color = term)) +
geom_point() +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
geom_vline(xintercept = 0)
o este:
mod_1 %>% augment(data = df_m) %>%
ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1.f)) +
geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
geom_line()
mod_1 %>% augment(data = df_m) %>%
ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1.f)) +
geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
geom_line() +
facet_wrap(vars(ESTUDIOM.f))
mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2, data = df_m)
mod_2 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)
anova(mod_1, mod_2)
#> Analysis of Variance Table
#>
#> Model 1: PESON1 ~ SEMANAS + EDADM.2
#> Model 2: PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS * EDADM.2)
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 57048 10755232331
#> 2 57047 10754104306 1 1128025 5.9838 0.01444 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(mod_1, mod_2)
#> df AIC
#> mod_1 4 854907.8
#> mod_2 5 854903.8
lmtest::lrtest(mod_1, mod_2)
#> Likelihood ratio test
#>
#> Model 1: PESON1 ~ SEMANAS + EDADM.2
#> Model 2: PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS * EDADM.2)
#> #Df LogLik Df Chisq Pr(>Chisq)
#> 1 4 -427450
#> 2 5 -427447 1 5.9839 0.01444 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(!!!) Vamos a ordenar los modelos en función de su AIC. Para ello vamos a crear una lista con modelos
modelos <- list(mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2, data = df_m),
mod_2 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m) )
modelos_ordered_AIC <- purrr::map_df(modelos, broom::glance, .id = "model") %>% arrange(AIC)
modelos_ordered_AIC %>% gt::gt()
| model | r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 0.2720188 | 0.2719805 | 434.1809 | 7105.43 | 0 | 4 | -427446.9 | 854903.8 | 854948.6 | 10754104306 | 57047 |
| 1 | 0.2719425 | 0.2719169 | 434.1999 | 10654.22 | 0 | 3 | -427449.9 | 854907.8 | 854943.6 | 10755232331 | 57048 |
broom (antes creamos una función (!!!!)):my_kable <- function(df){ gt::gt(mutate_if(df, is.numeric, round, 2)) }
tidy(mod_1) %>% my_kable
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -2327.62 | 39.76 | -58.55 | 0 |
| SEMANAS | 140.92 | 0.97 | 145.97 | 0 |
| EDADM.2 | 2.33 | 0.35 | 6.75 | 0 |
Con el paquete modelr podemos fácilmente comparar las predicciones de varios modelos:
zz <- df_m %>% modelr::gather_predictions(mod_1, mod_2)
Ejercicio: Pasad zz a formato ancho para ver las predicciones de los 2 modelos:
zz1 <- pivot_wider(zz, names_from = model, values_from = pred)
Por ejemplo un Logit:
mod_logit <- glm(CESAREA.f ~ PESON1 + SEMANAS + EDADM.2 + ESTUDIOM.f, family = binomial(link = "logit"), data = df_m)
summary(mod_logit)
#>
#> Call:
#> glm(formula = CESAREA.f ~ PESON1 + SEMANAS + EDADM.2 + ESTUDIOM.f,
#> family = binomial(link = "logit"), data = df_m)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -1.8803 -0.8662 -0.7599 1.3804 2.1197
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.23616479 0.20327735 11.001 < 0.0000000000000002 ***
#> PESON1 0.00013040 0.00002128 6.129 0.000000000882 ***
#> SEMANAS -0.13285957 0.00567450 -23.413 < 0.0000000000000002 ***
#> EDADM.2 0.04918850 0.00191326 25.709 < 0.0000000000000002 ***
#> ESTUDIOM.fMedios 0.10281855 0.02322645 4.427 0.000009564693 ***
#> ESTUDIOM.fUniversidad 0.07570366 0.02335789 3.241 0.00119 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 69663 on 57050 degrees of freedom
#> Residual deviance: 68142 on 57045 degrees of freedom
#> AIC: 68154
#>
#> Number of Fisher Scoring iterations: 4
En los modelos lineales, calcular efectos marginales es sencillo, pero el Logit es un modelo no lineal. Para calcular efectos marginales en modelos no lineales podemos usar el paquete margins. Por ejemplo, calculemos el Efecto marginal medio o average marginal effect (AME) en mod_logit.
mod_1_AME <- mod_logit %>% margins::margins() %>% summary()
mod_1_AME
#> factor AME SE z p lower upper
#> EDADM.2 0.0100 0.0004 26.1951 0.0000 0.0093 0.0108
#> ESTUDIOM.fMedios 0.0210 0.0047 4.4201 0.0000 0.0117 0.0303
#> ESTUDIOM.fUniversidad 0.0154 0.0047 3.2374 0.0012 0.0061 0.0247
#> PESON1 0.0000 0.0000 6.1353 0.0000 0.0000 0.0000
#> SEMANAS -0.0271 0.0011 -23.8118 0.0000 -0.0294 -0.0249
Si queremos calcular efectos marginales para unos valores concretos de los regresores
mod_logit %>% margins::margins(at = list(SEMANAS = c(25, 35), ESTUDIOM.f = c("Primarios", "Medios", "Universidad")), variables = "EDADM.2" )
Si prefieres visualizarlo, utiliza margins::cplot():
margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM.2", what = "effect", drop = TRUE)
margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM.2")
sjPLot y friends …m1 <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 + EDADP.2, data = df)
m2 <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 + EDADP.2 + ESTUDIOM.f + ESTUDIOP.f, data = df)
sjPlot::tab_model(m1)
sjPlot::plot_model(m1, sort.est = TRUE)
sjPlot::tab_model(m1, m2)
stargazermod_1 <- lm(PESON1 ~ SEMANAS + SEXO1.f , data = df_m)
stargazer::stargazer(mod_1, type = "html")
| Dependent variable: | |
| PESON1 | |
| SEMANAS | 141.401*** |
| (0.955) | |
| SEXO1.fMujer | -123.584*** |
| (3.604) | |
| Constant | -2,211.576*** |
| (37.258) | |
| Observations | 57,051 |
| R2 | 0.286 |
| Adjusted R2 | 0.286 |
| Residual Std. Error | 429.964 (df = 57048) |
| F Statistic | 11,430.040*** (df = 2; 57048) |
| Note: | p<0.1; p<0.05; p<0.01 |
modelsummary#remotes::install_github('vincentarelbundock/modelsummary')
library(modelsummary)
mys_modelitos <- list()
mys_modelitos[["PESON1: OLS 1"]] <- lm(PESON1 ~ SEMANAS + SEXO1.f, df_m)
mys_modelitos[["PESON1: OLS 2"]] <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 , df_m)
mys_modelitos[["CESAREA: Logit 1"]] <- glm(CESAREA.f ~ SEMANAS + SEXO1.f , data = df_m, family = binomial(link = "logit"))
mm <- msummary(mys_modelitos, title = "Resultados de estimación")
mm
| Resultados de estimación | |||
|---|---|---|---|
| PESON1: OLS 1 | PESON1: OLS 2 | CESAREA: Logit 1 | |
| (Intercept) | -2211.576 | -2301.307 | 3.761 |
| (37.258) | (39.360) | (0.184) | |
| SEMANAS | 141.401 | 141.711 | -0.118 |
| (0.955) | (0.956) | (0.005) | |
| SEXO1.fMujer | -123.584 | -123.749 | -0.061 |
| (3.604) | (3.602) | (0.018) | |
| EDADM.2 | 2.406 | ||
| (0.342) | |||
| Num.Obs. | 57051 | 57051 | 57051 |
| R2 | 0.286 | 0.287 | |
| Adj.R2 | 0.286 | 0.287 | |
| AIC | 853789.1 | 853741.5 | 69027.1 |
| BIC | 853824.9 | 853786.2 | 69054.0 |
| Log.Lik. | -426890.542 | -426865.740 | -34510.573 |
reportsreportslibrary(reports) #- devtools::install_github("neuropsychology/report")
my_model <- lm(PESON1 ~ SEMANAS + SEXO1.f, df_m)
rr <- report(my_model, target = 1)
rr
to_table(rr) %>% gt::gt()
to_fulltext(rr)
Recuerda también que se puede obtener la ecuación de un modelo con:
library(equatiomatic) #- remotes::install_github("datalorax/equatiomatic")
extract_eq(mod_1)
#> $$
#> \text{PESON1} = \alpha + \beta_{1}(\text{SEMANAS}) + \beta_{2}(\text{SEXO1.f}_{\text{Mujer}}) + \epsilon
#> $$
extract_eq(mod_1, use_coefs = TRUE)
#> $$
#> \text{PESON1} = -2211.58 + 141.4(\text{SEMANAS}) - 123.58(\text{SEXO1.f}_{\text{Mujer}}) + \epsilon
#> $$
Para que los resultados salgan rápido y se puedan visualizar, vamos a restringir el df:
df_x <- df %>% select(SEMANAS, CESAREA.f, EDADM.2, PESON1, SEXO1.f) %>% #- cogemos unas pocas variables
na.omit() %>% #- quitamos NA's
dplyr::sample_n(200) %>% #- y solo 200 bebes
mutate_if(is.factor, as.numeric) #- tienen q ser numericas
DataExplorer::plot_prcomp() para visualizar los resultados obtenidos con stats::prcomp()DataExplorer::plot_prcomp(df_x, variance_cap = 0.9, nrow = 2L, ncol = 2L)
FactoMineR::PCA()FactoMineR::PCA(df_x)
factoextra para visualizarres.pca <- FactoMineR::PCA(df_x, ncp = 5, graph = FALSE)
factoextra::fviz_screeplot(res.pca, addlabels = TRUE)
factoextra::fviz_pca_var(res.pca, col.var = "black")
factoextra::fviz_contrib(res.pca, choice = "var", axes = 1, top = 10) # Contributions of variables to PC1
factoextra::fviz_pca_ind(res.pca) #- resultados para bebes individuales
factoextra::fviz_pca_ind(res.pca, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)# Avoid text overlapping (slow if many points))
psycho::n_factors(): Find optimal components number using maximum method aggreement, pero no me acaba de funcionarlibrary(psycho)
results <- df_numeric %>% psycho::n_factors()
print(results)
plot(results)
summary(correlation(df_numeric))
df_xx <- base::scale(df_x)
k3 <- stats::kmeans(df_x, centers = 3, nstart = 7)
factoextra::fviz_cluster(k3, data = df_x)
Hice este gráfico para ver si lo entendía, pero …
df_x %>% as_tibble() %>% mutate(cluster = k3$cluster) %>%
ggplot(aes(SEMANAS, PESON1, shape = as.factor(CESAREA.f), color = factor(cluster), label = SEXO1.f)) +
geom_point() +
geom_text()
#- https://uc-r.github.io/kmeans_clustering
k2 <- kmeans(df_x, centers = 2, nstart = 7)
k3 <- kmeans(df_x, centers = 3, nstart = 7)
k4 <- kmeans(df_x, centers = 4, nstart = 7)
k5 <- kmeans(df_x, centers = 5, nstart = 7)
p1 <- fviz_cluster(k2, geom = "point", data = df_x) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = df_x) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = df_x) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = df_x) + ggtitle("k = 5")
library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow = 2)
#- https://uc-r.github.io/kmeans_clustering
distance <- factoextra::get_dist(df_x)
factoextra::fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
factoextra::fviz_nbclust(df_xx, kmeans, method = "gap_stat")
explore::explain_tree() que: Explain a target using a simple decision tree (classification or regression)df_m %>% explore::explain_tree(target = PESON1)
df_m %>% explore::explain_tree(target = CESAREA.f)
df_m %>% explore::explain_tree(target = SEXO1.f)
df_m %>% explore::explain_tree(target = ESTUDIOM.f)
fit <- rpart::rpart(SEXO1.f ~ . , data = df_m)
rattle::fancyRpartPlot(fit, palettes = c("Greens", "Reds"), sub = "")
Simplemente algunas referencias:
Machine Learning for Everyone. Un post introductorio.
Hands-On Machine Learning with R. Un buen bookdown
101 Machine Learning Algorithms for Data Science with Cheat Sheets. Cheatsheet de algoritmos ML.
101 Data Science Interview Questions, Answers, and Key Concepts. Post con cosas que “deberías” saber si quieres un trabajo como Data Scientist.
Machine learning essentials. Un tema de un curso sobre ML.
Introduction to Machine Learning with R. Otro curso sobre ML.
Machine Learning. Bookdown centrado en explicar las principales ideas/conceptos de ML.
Machine Learning con R y caret. Una buena introducción A R con caret. En castellano, por Joaquín Amat Rodrigo. En su página web tiene muchos recursos interesantes sobre estadística y ML.
Customer Churn Modeling using Machine Learning with parsnip. Un post con un buen ejemplo de uso de ML.
Choosing the right estimator. Guía de scitit learn para selección de modelo.
Una introducción visual al machine learning - I. Post introductorio pero muy bonito visualmente sobre ML. Aquí la segunda parte.
How to use recipes package from tidymodels for one hot encoding
R: MLR, Decision Trees and Random Forest to Predict MPG for 2019 Vehicles. Después se fueron a hacerlo en TensorFlow
El warning nos salta porque algunas variables son factores y tienen NA’s, pero NA no es uno de los levels del factor. El warning simplemente nos recuerda que hay como un level “oculto”.↩
Por defecto, el paquete usa Eicker-Huber-White robust standard errors, habitualmente conocidos como errores estándar “HC2”. Se pueden especificar otros métodos con la opción se_type. Por ejemplo se puede utilizar el método usado por defecto en Stata. Aquí puedes encontrar porque los errores estándar de Stata difieren de los usados en R y Phyton↩